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The  elastic  constants  of  a  wide  range  of  models  of  defected  crystalline  and  amorphous  silicon  are  calculated, 
using  the  environment-dependent  interatomic  potential  (EDIP).  The  defected  crystalline  simulation  cells  con¬ 
tain  randomly  generated  defect  distributions.  An  extensive  characterization  of  point  defects  is  performed, 
including  structure,  energy  and  influence  on  elastic  constants.  Three  important  conclusions  are  drawn.  (1) 
Defects  have  independent  effects  on  the  elastic  constants  of  silicon  up  to  (at  least)  a  defect  concentration  of 
0.3%.  (2)  The  linear  effect  of  Frenkel  pairs  on  the  (110)  Young’s  modulus  of  silicon  is  -1653  GPa  per  defect 
fraction.  (3)  17  different  point  defect  types  cause  a  very  similar  decrease  in  the  (110)  Young's  modulus: 
-(0.28  +  0.05)%  when  calculated  in  isolation  using  a  1728-atom  cell.  These  principles  will  be  very  useful  for 
predicting  the  effect  of  radiation  damage  on  the  elastic  modulus  of  silicon  in  the  typical  case  in  which 
point-defect  concentrations  can  be  estimated,  but  the  exact  distribution  and  species  of  defects  is  unknown.  We 
also  study  amorphous  samples  generated  in  quenching  the  liquid  with  EDIP,  including  an  ideal  structure  of 
perfect  fourfold  coordination,  samples  with  threefold  and  fivefold  coordinated  defects,  one  with  a  nanovoid, 
and  one  with  an  amorphous  inclusion  in  a  crystalline  matrix.  In  the  last  case,  a  useful  Ending  is  that  the  change 
in  the  Young’s  modulus  is  simply  related  to  the  volume  fraction  of  amorphous  material,  as  has  also  been 
observed  by  experiment. 
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I.  INTRODUCTION 

Defects  in  silicon  have  been  more  extensively  studied  for 
their  effects  on  electronic  properties  than  for  their  effects  on 
mechanical  properties.  This  is  not  surprising,  given  the  ex¬ 
tensive  uses  that  have  been  made  of  silicon’s  electronic  prop¬ 
erties.  With  the  advent  of  ever  more-sensitive  microelectro¬ 
mechanical  (MEMS)  devices,  however,  the  importance  of 
precise  knowledge  of  the  mechanical  properties  of  compo¬ 
nent  materials  has  grown.  Shifts  in  mechanical  properties 
may  well  compromise  the  functioning  of  a  highly  sensitive 
MEMS  device.  Radiation  damage  and  even  small  changes  in 
temperature  or  stress  state  can  cause  sufficient  alteration  of 
dimension  and  elasticity  to  be  of  concern. 

Radiation-induced  changes  in  the  mechanical  properties 
of  silicon,  a  common  MEMS  material,  have  been  studied  for 
rather  large  ion  fluences  and  elastic  constant  changes,  from 
S5%  (Ref.  1)  to  complete  amorphization.2  4  Such  high  lev¬ 
els  of  radiation  damage  are  common  when  performing  ion 
implantation.  It  is  not  clear,  however,  that  effects  at  very  low 
radiation  doses,  such  as  might  occur  in  an  environment  such 
as  space,  can  be  correctly  extrapolated  from  such  high-dose 
regimes.  While  not  simulating  radiation  damage  per  se,  the 
present  study  addresses  elastic  constant  changes  occurring  in 
defected  and  amorphous  silicon,  both  possible  results  of  ir¬ 
radiation. 

Molecular  dynamics  simulations  predict  that  radiation 
damage  consists  of  a  mixture  of  isolated  defects,  aggregate 
defects,  and  amorphous  regions. 5-8  Whether  the  isolated  de¬ 
fects  or  amorphous  regions  dominate  will  depend  on  such 


factors  as  the  incident  particle  flux,  mass,  and  energy,  as  well 
as  the  temperature  of  the  damaged  material.  Heavier,  more 
energetic  recoils  produce  larger  amorphous  pockets  which 
are  unlikely  to  anneal,  even  at  room  temperature.5  The  quali¬ 
tative  results  of  these  studies  are  relatively  insensitive  to  the 
interatomic  potentials  used,  among  the  most  popular  and 
well-tested  models  of  silicon:  the  Stillinger- Weber  (SW) 
potential9  (Refs.  5,  7,  and  8),  the  Tersoff  potentials10,11  (Refs. 
6-8),  and  the  environment-dependent  interatomic  potential 
(EDIP)12,13  (Ref.  8). 

For  the  present  study  of  crystalline  defects  and  amorphous 
silicon,  EDIP  is  a  natural  choice.  It  was  fit  to  a  few  point 
defect  energies  and  then  tested  for  a  wide  variety  of  unre¬ 
lated  structures,  with  considerable  success  in  light  of  its  rela¬ 
tively  small  number  of  parameters  (13).  For  example,  among 
the  most  commonly  used  potentials,  EDIP  is  the  only  one  to 
correctly  predict  dislocation  core  reconstructions  and  a  direct 
quench  from  the  liquid  to  a  high  quality  amorphous  phase.13 
Since  the  original  study,  EDIP  has  been  used  extensively  in 
simulations  of  crystalline  defects14-19  and  amorphous 
structures,20-25  so  it  may  be  considered  well  tested  for  the 
present  application. 

More  importantly,  from  a  physical  point  of  view,  EDIP 
includes  environment-dependent  changes  in  chemical  bond¬ 
ing,  inferred  directly  from  ab  initio  calculations  and  experi¬ 
mental  data. 12,26,27  Specifically,  the  bond  order  (strength  of 
the  pairwise  attraction),  the  preferred  bond  angle,  and  the 
angular  stiffness  depend  strongly  on  the  local  coordination 
number  Z.  In  contrast,  the  form  of  the  SW  potential  can  only 
be  justified  for  rigid  sp3  hybrid  bonds,  which,  in  reality,  are 
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replaced  by  other  covalent  hybrids  or  metallic  bonds  in  non- 
tetrahedral  defects  and  disordered  structures.12,26  The  reason¬ 
able  description  of  the  metallic  liquid  and  certain  defects  by 
SW,  therefore,  should  be  viewed  as  fortuitous.  In  compari¬ 
son,  the  Tersoff  potentials  incorporate  more  realistic  features 
of  the  bond  order,  which  can  be  derived  analytically  from 
tight-binding  models.28  The  functional  form,  however,  is  in¬ 
consistent  with  silicon  elastic  constant  relations12  (satisfied 
by  SW  and  EDIP)  and  seems  unable  to  simulataneously  de¬ 
scribe  elasticity,  defects,  and  phase  transitions.29 

The  importance  of  the  choice  of  potential  is  well  illus¬ 
trated  by  the  present  topic.  Experiments  and  simulations  both 
show  that  amorphous  silicon  is  less  stiff  than  the 
crystal.2,3,30,31  On  the  other  hand,  the  only  prior  work  on  the 
effect  of  defects  at  low  concentrations  on  crystal  elasticity  by 
Clark  and  Ackland  (CA)  concludes  that  both  vacancies  and 
interstitials  tend  to  increase  the  elastic  constants,32  although 
one  might  expect  the  opposite  trend,  since  the  accumulated 
effect  of  defects  in  the  crystal  roughly  approaches  an  amor¬ 
phous  structure.  These  authors  introduced  a  new  potential 
depending  only  upon  pairwise  bond  lengths,  which  has  since 
received  very  little  testing.  A  serious  concern  regarding  their 
results  is  the  poor  description  of  the  crystal  elastic  constants 
by  the  CA  potential  (compared  to  SW  and  EDIP),  which  can 
be  attributed  to  the  lack  of  explicit  angular  dependence,  in 
violation  of  silicon  elastic  constant  relations12  and  analysis  of 
tight-binding  models.28,33  Following  Ackland,34,35  the  CA  po¬ 
tential  also  assumes  exactly  four  bonds  per  atom,  which  ob¬ 
viously  cannot  apply  to  most  defects.  Therefore,  it  is  not 
surprising  that  in  the  present  study  with  EDIP  we  reach  a 
very  different  conclusion:  point  defects  tend  to  make  the 
crystal  more  soft  (as  is  the  amorphous  phase).  We  check  that 
the  same  trend  is  predicted  by  SW,  so  it  is  clear  that  angular 
terms,  neglected  by  CA,  play  an  important  role. 

In  practical  situations,  it  is  the  finite-temperature  behavior 
of  silicon  elasticity  that  would  be  of  interest.  While  statistical 
mechanical  methods  do  exist  for  calculating  finite- 
temperature  elastic  constants  (e.g.,  Monte  Carlo36  or  molecu¬ 
lar  dynamics37  simulations),  we  decided  not  to  pursue  such 
methods  in  this  study.  We  justify  this  choice  in  light  of  the 
small  defect  fractions  and  concurrently  small  shifts  in  elastic 
constant  that  were  considered  here,  which  would  have  been 
difficult  to  notice  with  the  slower  convergence  of  fluctuation 
methods.  In  addition  to  having  greater  precision,  the  cell 
deformations  described  below  were  much  simpler  and  com¬ 
putationally  cheaper  than  finite-temperature  methods.  Since 
the  elastic  constants  of  silicon  shift  only  by  around  1%  be¬ 
tween  0  and  300  K,38  we  believe  that  the  trends  in  the  elastic 
constants  with  defect  content  are  essentially  captured  by  our 
zero-temperature  calculations. 

II.  METHODS 
A.  Sample  creation 
1.  Defected  samples 

Defected  samples  were  created  by  inserting  vacancies  and 
interstitials  into  1728-atom  supercells  (6X6X6  unit  cells), 
using  periodic  boundary  conditions.  Randomly  selected  at¬ 


oms  were  removed  to  produce  vacancies  and  random  posi¬ 
tions  were  chosen  to  insert  extra  atoms.  Three  types  of 
sample  were  produced:  one  containing  only  vacancies,  one 
having  only  interstitials,  and  one  containing  Frenkel  pairs. 
Fifteen  samples  were  produced  of  each  sample  type,  consist¬ 
ing  of  three  sets  of  supercells  containing  1,2, 3,4,  and  5  de¬ 
fects,  for  a  grand  total  of  45  sample  supercells. 

Upon  insertion  of  defects,  the  supercells  were  relaxed  at 
0  K  by  relaxing  atomic  positions  while  iterating  over  cell 
dimension  changes.  Isotropic  expansion  and  cell-length 
changes  along  individual  axes  were  iteratively  explored  fol¬ 
lowing  a  simple  energy  gradient  algorithm.  The  supercells 
were  then  annealed  at  300  K  for  500  ps  using  the  DL  POLY 
(Ref.  39)  NPT  ensemble  with  Berendsen  thermostat  and 
barostat.  After  annealing,  the  supercells  were  again  relaxed  at 
0  K. 

It  was  hoped  that  this  procedure  would  yield  samples  rep¬ 
resenting  several  different  defect  types  in  random  geom¬ 
etries,  as  might  result  in  regions  of  a  collision  cascade  rich  in 
point  defects.  Some  samples  were  thrown  out  and  regener¬ 
ated,  since  sometimes  a  randomly  chosen  vacancy  and  a  ran¬ 
domly  placed  extra  atom  annihlated  either  immediately  or 
upon  annealing.  Other  times,  a  randomly  chosen  position  for 
an  extra  atom  was  too  close  to  a  crytalline  position  in  the 
supercell,  and  mayhem  resulted  upon  relaxation  due  to  the 
large  forces  present.  Samples  were  visually  inspected  to  en¬ 
sure  that  the  desired  scattered  defects  were  present  and  dis- 
cernably  distinct.  In  most  cases,  only  isolated  defects  were 
present.  Exceptions  included  a  two-interstitial  agglomerate 
( H 2) ,  and  divacancy  complexes  involving  missing  nearest 
neighbors  (V^/v)  and  missing  next-nearest  neighbors  (V^aw)- 
A  breaking  of  symmetry  in  the  (IAaw)  complex  yielded  an¬ 
other  divacancy  type  (V^J- 

In  addition  to  these  45  defected  samples,  17  additional 
1728-atom  supercells  were  prepared  in  order  to  characterize 
the  formation  energy  (£)),  volume  (Vj),  and  effects  on  elastic 
constants  of  various  interstitial  and  vacancy  configurations. 
These  samples  were  created  by  intentionally  arranging 
atomic  positions  within  the  supercell  then  relaxing  at  0  K,  as 
above.  Ej,  Vf,  and  elastic  constants  for  each  sample  were 
calculated  before  annealing  the  sample  at  300  K  for  500  ps, 
as  above.  Calculations  were  then  repeated  on  the  annealed 
samples. 

2.  Amorphous  samples 

Brief  descriptions  of  the  defect  content  of  amorphous 
samples  are  given  in  Table  I.  Amorphous  samples  A,  D,  and 
G  were  prepared  by  quenching  from  the  liquid  at  zero  pres¬ 
sure  as  follows.  A  diamond  crystal  was  melted  at  3000  K  for 
50  ps,  cooled  to  1500  K  over  100  ps,  then  equilibrated  at 
1500  K  for  an  additional  100  ps.  The  samples  were  then 
cooled  1000  K  over  1  ns  (the  transition  from  liquid  to  the 
amorphous  state  occurred  at  this  step),  then  annealed  at 
1000  K  for  2  ns.  The  annealed  amorphous  sample  was  then 
cooled  to  0  K  over  2  ns. 

Amorphous  samples  B  and  C  were  derived  from  an  inter¬ 
mediate  sample  that  was  prepared  from  sample  A  as  follows. 
A  negative  pressure  of  -100  GPa  was  imposed  on  sample  A 
at  0  K,  then  the  sample  was  annealed  at  1000  K  for  2  ns  and 


134113-2 


PHYSICAL  REVIEW  B  70,  134113  (2004) 


ELASTIC  CONSTANTS  OF  DEFECTED  AND  AMORPHOUS... 


TABLE  I.  Descriptions  of  amorphous  samples.  Calculation  of 
numbers  of  fivefold  and  threefold  coordination  defects  present  used 
a  rounding  of  the  EDIP  coordination  value  Z. 


Sample 

Size 

Fivefold 

Threefold 

Other 

A 

216 

6 

0 

B 

216 

5 

1 

C 

216 

13 

3 

Void  present 

D 

216 

10 

0 

E 

1728 

116 

2 

F 

1728 

112 

2 

G 

64 

0 

0 

slowly  cooled  to  0  K  at  constant  volume.  Finally,  the  struc¬ 
ture  was  relaxed  to  zero  pressure. 

Sample  B  was  derived  from  this  intermediate  structure  by 
annealing  at  1 100  K  for  2  ns  and  cooling  to  0  K,  all  at  zero 
pressure.  Sample  C  was  derived  from  the  same  intermediate 
structure  by  annealing  at  1 100  K  for  4  ns  and  then  cooling  to 
0  K,  all  at  constant  volume.  The  defect  content  of  sample  B 
was  similar  to  that  of  sample  A,  whereas  sample  C  contained 
a  sizable  void. 

Two  1728-atom  amorphous  samples  were  prepared  by 
quenching  from  high  temperature  at  constant  pressure. 
Sample  E  was  prepared  by  heating  a  crystalline  sample  to 
3000  K,  equilibrating  for  300  ps,  then  cooling  to  10  K  at 
1  K/ps.  Sample  F  started  with  random  atomic  positions  at 
5000  K,  then  was  cooled  to  10  K,  also  at  1  K/ps.  All  amor¬ 
phous  samples  used  here  were  annealed  at  300  K  for  500  ps, 
then  relaxed  at  0  K  as  for  the  defected  samples  above. 

3.  Amorphous  pocket  sample 

A  composite  sample  consisting  of  an  amorphous  block 
surrounded  by  crytalline  material  was  prepared  by  embed¬ 
ding  sample  G  into  a  crystalline  matrix.  This  was  done  by 
first  cutting  out  a  2X2X2  unit  cell  cube  from  a  6X6X6 
unit  cell  supercell.  The  crystalline  supercell  was  then  homo¬ 


geneously  and  anisotropically  expanded  to  accommodate  the 
relaxed  sample  G.  The  sample  containing  the  amorphous 
pocket  was  then  relaxed  at  0  K  as  above,  then  annealed  as 
above  at  300  K  and  relaxed  again  at  OK. 


B.  Elastic  constant  calculation 
1.  Approach 

Strain  can  be  defined  as40 

eij=\(eij  +  ejl),  (1) 

where 


dllj 

eu  =  — 
1  dX: 


(2) 


with  u  being  displacement  and  x  being  position.  The  work 
needed  to  impose  a  strain  e  is  (switching  to  Voigt  notation) 


A  E 
V 


2C‘J£ieP 


(3) 


where  C(/-  are  the  elastic  constants  and  V  is  the  volume.  For  a 
cubic  system  such  as  silicon,  only  three  independent  con¬ 
stants  exist.  We  make  the  approximation  for  our  defected 
samples  that  this  symmetry  remains  intact. 

Computationally  convenient  strains  to  impose  on  an 
orthorhombic  supercell  aligned  with  the  diamond  unit  cell 
axes  are  suggested  by  Eq.  (3).  The  tensor  strain  notation 
shown  in  Eq.  (4)  is  used  in  Eqs.  (5)-(7),  with  nonspecified 
strains  all  equal  to  zero: 


1 

1 

G  2  e6 

2f5 

1 

1 

V6 

6 A 

2  4 

1  1 

€a 

2  5  2  4 

e3 

TABLE  II.  Analytical  and  numerically  (energy  curve  fitting)  computed  EDIP  elastic  constants  (GPa)  and 
Kleinman’s  internal  strain  parameter  £  for  the  diamond  phase.  C44  is  the  shear  modulus  when  internal 
relaxation  is  disallowed.  The  values  published  by  Justo  et  al.  are  provided  here  for  comparison. 


Analytical 

Numerical 

Justo  et  al.  (Ref.  13) 

Expt.  (Refs.  38  and  41) 

Cu 

171.99 

172.00 

175 

167.54 

Cl2 

64.72 

64.71 

62 

64.92 

C44 

72.75 

72.75 

71 

80.24 

c° 

^44 

112.39 

112.39 

112 

IIId 

B 

100.47 

100.47 

99 

99.13 

i 

0.51727 

0.51727“ 

0.497b 

0.72 

£<100>C 

136.60 

136.62 

143 

131.28 

£<no>c 

164.05 

164.05 

164 

170.63 

“Calculated  computationally  by  shearing  a  supercell  and  measuring  atomic  relaxation. 
bCalculated  from  C,j  (Ref.  27). 

Calculated  from  CtJ  (Ref.  40). 

dNot  experimentally  accessible.  Calculated  using  ab  initio  (LDA)  methods  in  Ref.  12. 
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Cn  = 


1  &E 
Vd^’ 


C\\  +  C 12- 


1  &E 

^  —  ^1  —  ^9’ 

2V  de1 


(5) 

(6) 


_  1  d1 E 

°44=v^- 


(7) 


By  imposing  isotropic  strains  (el  =  e2=e3),  we  can  also  ob¬ 
tain  the  bulk  modulus 


B  =  V 


_f?E 

dV2' 


(8) 


2.  Defected  samples 

Elastic  constants  were  calculated  by  fitting  polynomials  to 
plots  of  energy  vs  strain.  Excellent  results  were  obtained  for 
a  wide  range  of  energy  sampling  intervals  in  strain  for  the 
crystal  and  for  samples  containing  only  point  defects,  but  the 
amorphous  samples  were  found  to  be  apt  to  change  bonding 
structure  with  large  imposed  strains.  By  trial  and  error,  it  was 
determined  that  a  step  size  in  strain  of  1  X  10-5  was  suffi¬ 
ciently  small  to  not  cause  bonding  rearrangements  in  the 
amorphous  samples,  yet  produced  large  enough  energy  dif¬ 
ferences  to  numerically  extract  second  derivatives  from  fit 
polynomials. 

The  number  of  points  to  fit  was  also  decided  by  experi¬ 
mentation.  The  second  deriviative  of  the  energy  was  of  pri¬ 
mary  concern  here.  A  greater  number  of  points  would  help  to 
recognize  and  account  for  higher-order  effects,  but  a  large 
number  of  points  at  large  strains  would  also  allow  higher- 
order  effects  to  dominate  the  fit.  A  total  of  21  points  in  the 
strain  interval  +1  X  1  (L 4  was  calculated  for  each  sample  dis¬ 
cussed  here,  but  it  was  found  that  fitting  only  the  nine  middle 
points  yielded  the  most  consistent  results. 

Starting  from  the  relaxed  structure,  strains  were  sequen¬ 
tially  imposed  in  steps  of  1  X  10-5.  Atomic  positions  were 
relaxed  at  each  strain  step.  Starting  again  from  zero  strain, 
negative  strain  was  then  incrementally  imposed.  Fourth  and 
fifth  order  fits  to  the  resultant  curves  were  found  to  give 
essentially  identical  results,  with  lower-order  fits  giving 
slightly  different  results,  and  higher-order  fits  being  numeri¬ 
cally  unreliable.  All  results  reported  here  are  based  on  fourth- 
order  least-squares  fits. 

To  help  minimize  the  error  in  our  assumption  that  cubic 
symmetry  holds  for  the  defected  samples,  all  C,;  were  inde¬ 
pendently  calculated  for  all  three  orthogonal  directions,  then 
averaged.  The  variation  of  C44  among  the  three  directions 
was  greatest,  and  rose  to  ~  1  %  for  a  few  of  the  samples,  but 
was  typically  similar  to  that  of  Cn  and  C12,  namely,  ;£0.1%. 
The  lesser  precision  of  C44  is  to  be  expected,  given  the 
sample  preparation  method  described  in  Sec.  II  A  1,  which 
guarantees  that  the  relaxed  sample  will  occupy  an  energy 
minimum  with  respect  to  the  isotropic  and  monoaxial  expan¬ 
sions  used  to  obtain  energy  curves  for  B  and  Cn,  but  will  not 
guarantee  this  for  the  shear  applied  to  calculate  C44. 


0.0  0.1  0.2  0.3 


Percent  Defect 


FIG.  1.  Elastic  constants  as  a  function  of  point  defect  fraction. 
Diamonds  indicate  interstitials,  triangles  vacancies,  and  squares 
Frenkel  pairs.  Each  point  represents  the  average  of  three  samples. 
Bars  indicate  the  standard  deviation  of  the  three  samples.  Lines  are 
unweighted  least  squares  fits  constrained  to  match  the  crystalline 
value  on  the  Y  axis. 


As  a  validation  of  our  numerical  approach,  we  compared 
values  found  for  the  crystal  with  those  previously  published 
for  EDIP.13  When  discrepancies  were  observed,  we  calcu¬ 
lated  analytical  values27  for  the  EDIP  diamond  lattice  elastic 
constants,  and  found  that  our  numerically  calculated  values 
agreed  to  within  about  four  significant  figures.  The  results  for 
these  values  are  shown  in  Table  II.  It  can  be  seen  that  the 
values  of  Justo  et  al.,  calculated  using  a  deformation/energy 
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TABLE  III.  Fit  slopes  to  the  points  in  Figs.  1  and  2.  Units  are  per  defect  fraction,  with  elastic  constants  in  GPa.  “Difference”  refers  to 
that  between  the  Frenkel  pair  result  and  the  sum  of  the  vacancy  and  interstitial  results.  The  numbers  in  parentheses  are  from  fits  to  the  data 
in  Clark  and  Ackland  (Ref.  32). 


Cn 

Cl2 

C44 

B 

-£(110) 

AVIV 

Vacancies 

-1030(821) 

-381(565) 

-369(-8.7) 

-593(648) 

-899(217) 

0.42 

Interstitials 

-89.2(715) 

733(443) 

— 440( — 57.9) 

462(531) 

-751(161) 

1.35 

Frenkel  pairs 

-1133 

339 

-809 

-145 

-1653 

1.78 

Vac.  +  Int. 

-1119 

352 

-810 

-130 

-1651 

1.77 

%  Difference 

-1.2 

3.6 

0.1 

-10.3 

-0.2 

-0.9 

method,  very  similar  to  our  results  in  the  “numerical”  col¬ 
umn  of  Table  H,  were  not  accurate  to  the  number  of  places 
published.  They  should  be  disregarded.  Details  on  calculat¬ 
ing  the  analytical  elastic  constants  from  the  potential  are 
given  in  Ref.  27. 


3.  Amorphous  samples 

In  the  case  of  a  truly  amorphous  material,  the  elastic  prop¬ 
erties  will  be  isotropic,  and  only  two  independent  elastic 
constants  will  exist.  The  finite  size  and  periodic  boundary 
conditions  imposed  on  the  samples  of  our  study  prohibit  this 
isotropic  ideal.  We  nevertheless  calculated  the  C-,j  for  our 
amorphous  Si  samples  as  above,  averaging  over  results  from 
three  independent  directions. 

An  additional  step  we  applied  to  the  amorphous  sample 
elastic  constant  results  was  to  distill  the  set  of  four  elastic 
constants  to  a  self-consistent  set  by  means  of  an  iterative 
process.  The  bulk  modulus,  being  direction-independent,  was 
left  alone.  By  pairing  each  of  the  three  C,j  with  B,  a  predic¬ 
tion  of  the  other  two  Cjf  s  was  made.  At  each  iteration,  the 
original  Ci;-  was  averaged  with  the  two  predicted  values  until 
a  self-consistent  set  was  achieved.  For  the  largest  samples, 
the  change  in  the  Cl;’s  due  to  this  process  was  sSO.  I  %  for  Cn 
and  C12,  and  about  0.5%  for  C44.  For  the  smallest  samples, 
these  changes  were  ~1  and  ~10%,  respectively.  The  small 
adjustments  necessary  to  achieve  self-consistency  of  the 
large  (1728  atoms)  samples  demonstrates  their  close  approxi¬ 
mation  to  the  isotropic  ideal. 

4.  Amorphous  pocket  sample 

One  of  the  goals  of  our  work  was  to  see  what  effect  amor¬ 
phous  regions  embedded  in  a  crystalline  matrix  (as  might  be 
residual  from  a  collision  cascade)  would  have  on  the  overall 
elastic  modulus.  As  for  the  defected  samples,  cubic  symme¬ 
try  doesn’t  strictly  apply  for  the  amorphous  pocket  sample 
studied  here,  nor  is  the  sample  isotropic.  We  wish  neverthe¬ 
less  to  compare  the  elastic  modulus  of  the  composite  with 
that  of  its  component  parts  (crystalline  and  amorphous).  We 
therefore  calculated  the  Cy’s  of  the  composite  as  for  the  de¬ 
fected  samples.  We  then  reduced  the  crystalline  and  compos¬ 
ite  Cjj  s  to  an  effective  isotropic  Young’s  modulus  E  by  using 
the  mean  of  the  Reuss  and  Voigt  spatial-average  limits.42 

For  isotropic  materials,  upper  and  lower  bounds  on  the 
Young’s  modulus  of  a  two-phase  composite  ( Ec )  are  given 
by4-43 


Vi  (1  -  Vi) 

£j+ E2 


<Ec<V1El  +  (l-Vl)E2, 


(9) 


where  V1  and  /:)  are,  respectively,  the  volume  fraction  and 
Young’s  modulus  of  the  embedded  material,  and  E2  is  the 
Young’s  modulus  of  the  matrix.  Clearly,  the  crystalline  ma¬ 
trix  we  had  here  is  not  isotropic,  but  the  Reuss  and  Voigt 
limits  for  polycrystalline  aggregates  provide  a  convenient 
way  to  test  the  applicability  of  Eq.  (9)  to  our  embedded 
amorphous  Si  sample. 


III.  RESULTS 
A.  Defected  samples 

The  results  of  the  elastic  constant  calculations  for  the 
samples  generated  by  random  placement  of  vacancies  and 
interstitials  are  summarized  in  Fig.  1 .  It  can  be  seen  that  the 
trend  of  the  elastic  constants  with  defect  content  is  approxi¬ 
mately  linear.  An  unweighted  fit  was  used,  since  the  small 
standard  deviation  of  the  single  defect  samples,  as  well  as  the 
small  number  of  samples,  made  weighted  fitting  difficult  to 
interpret.  The  fit  slopes  are  summarized  in  Table  III. 

The  isolated  defects  appear  to  be  having  largely  indepen¬ 
dent  effects  on  the  elastic  constants  and  volume,  as  evi¬ 
denced  by  the  linear  trends  in  Figs.  1  and  2,  as  well  as  by  the 
close  agreement  between  the  slope  of  the  Frenkel  pair  plots 
and  the  sum  of  the  vacancy  and  interstitial  slopes  as  shown 
in  Table  III.  Both  interstitials  and  vacancies  are  shown  to 
cause  volume  expansion  in  Fig.  2,  vacancies  to  a  lesser  de¬ 
gree  than  interstitials.  The  general  softening  of  the  elastic 
constants  with  rising  defect  concentration  (with  the  excep¬ 
tion  of  interstitials  for  C12  and  the  bulk  modulus,  and  Frenkel 


FIG.  2.  Percent  volume  change  as  a  function  of  point  defect 
fraction.  Symbols  are  as  in  Fig.  1. 
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TABLE  IV.  Elastic  constants  (GPa)  of  EDIP-generated  amorphous  samples  described  in  the  text.  Also  shown  are  percent  volume  change 
and  energy  gain  per  atom  with  respect  to  the  crystalline  values  (20.018  A3  atomic  volume  and  -4.6500  eV).  Results  from  experiment  (room 
temperature),  tight-binding  molecular  dynamics  (room  temperature),  and  EDIP  are  shown  for  comparison.  Where  published  results  were  not 
completely  consistent  with  isotropicity,  the  range  of  possible  derived  values  is  given. 


Cn 

^12 

C44 

B 

E 

%  AV 

A E,  eV 

A 

131.3 

80.4 

25.5 

91 A 

70.3 

3.8 

0.1885 

B 

130.2 

82.9 

23.6 

98.7 

65.7 

4.4 

0.1999 

C  (void) 

98.7 

57.9 

20.4 

71.5 

55.9 

10.9 

0.2432 

D 

134.7 

82.7 

26.0 

100.0 

71.7 

3.3 

0.1868 

E 

133.0 

81.9 

25.5 

99.0 

70.5 

3.1 

0.2059 

F 

132.9 

81.1 

25.9 

98.4 

71.4 

3.3 

0.2064 

G  (defect  free) 

131.0 

81.4 

24.8 

97.9 

68.6 

5.2 

0.1815 

EDIP  (Ref.  49) 

130 

81 

56 

97 

67-145 

2 

0.199 

EDIP  quench  (Ref.  8) 

3.5 

0.25 

EDIP  irradiated  (Ref.  8) 

2. 2-3. 6 

0.19-0.40 

Expt.  (Ref.  45)a 

156 

58.4 

48.8 

90.9 

124 

Expt.  (Ref.  2) 

156 

57.8 

49.2 

90.6 

125 

1.3 

Expt.  (Ref.  3) 

138 

42 

48 

74 

118 

TBMD  (Ref.  45) 

149 

46.9 

55.4 

75-84 

127-136 

Expt.  (Ref.  44) 

1.8 

“The  values  for  C/;  here  are  based  on  the  measured  Young’s  modulus  of  Tan  et  al.  (Ref.  31)  combined  with  a  Raleigh  wave  measurement 
discussed  in  De  Sandre  et  al.  (Refs.  45-48). 


pairs  for  C12)  does  not  agree  with  the  previous  finding  that 
both  interstitials  and  vacancies  in  silicon  stiffen  the  crystal. 
Vacancies  are  found  here  to  soften  each  of  the  elastic  con¬ 
stants  considered,  the  opposite  being  true  for  all  but  C44  ac¬ 
cording  to  Clark  and  Ackland.32  Points  of  agreement  be¬ 
tween  that  study  and  the  present  one  include  the  softening  of 
C44  by  both  vacancies  and  interstitials  (though  the  effect  is 
an  order  of  magnitude  greater  here),  and  the  stiffening  of 
both  B  and  C12  by  interstitials. 

B.  Amorphous  samples 

Elastic  constants  of  our  amorphous  samples  are  shown  in 
Table  IV,  along  with  experimental2,3,3 1,44-48  and 
calculated8,45,49  data  for  comparison.  It  can  be  seen  that  the 
differences  among  samples  is  slight,  with  the  exception  of 
sample  C,  which  contained  a  void.  The  elastic  constants  for 
our  amorphous  samples  appear  to  be  fairly  insensitive  to  the 
coordination  defects  described  in  Table  I.  None  of  elastic 
constants  of  sample  G,  which  has  perfect  four-coordination 
throughout,  is  an  outlier  when  compared  to  the  other 
samples.  Not  surprisingly,  sample  G  is  the  lowest  energy 
configuration.  It  is  interesting  to  note,  however,  that  the  vol¬ 
ume  change  of  sample  G  is  the  greatest  of  all  EDIP  samples 
reported  in  Table  IV  (with  the  exception,  of  course,  of 
sample  C). 

We  note  that  depsite  close  agreement  of  our  large  amor¬ 
phous  samples  with  the  amorphous  sample  of  Vink  et  al.49 
for  the  values  of  Cn  and  C12,  our  value  for  C44  is  about  a 
factor  of  2  different.  As  mentioned  above,  our  large  (1728 
atoms)  amorphous  samples  came  very  close  to  satisfying  the 
isotropic  ideal  represented  by 


C44=!(Cu-C12).  (10) 

Since  even  our  smallest  (64  atoms)  samples  came  within 
10%  of  satisfying  Eq.  (10),  with  all  larger  samples  being 
even  closer  than  this  (even  sample  C  with  its  void),  one 
would  expect  the  1000  atom  sample  of  Vink  et  al.  to  exhibit 
fairly  isotropic  properties.  Clearly,  the  result  of  Vink  et  al. 
for  C44  cannot  be  correct,  since  it  does  not  agree  with  Cu 
and  C12  according  to  Eq.  (10).  We  suggest  that  a  factor  of 
two  error  might  have  been  introduced  when  calculating  C44, 
which  would  mean  that  the  actual  C44  of  the  Vink  et  al. 
amorphous  sample  should  be  28  GPa,  closer  to  the  32  GPa 
that  Vink  et  al.  report  for  the  amorphous  silicon  modeled 
using  the  Stillinger- Weber  potential.49 

Although  the  bulk  modulus  is  correctly  predicted  to  be 
very  close  to  that  of  the  crystal,  the  overall  description  of 
amorphous  elasticity  is  not  completely  satisfactory.  Despite 
the  fact  that  amorphous  samples,  both  real  and  computa¬ 
tional,  have  properties  that  vary  according  to  the  method  of 
preparation,8  it  seems  likely  from  the  results  in  Table  IV  that 
EDIP  errs  systematically  in  predicting  the  elastic  constants  of 
amorphous  silicon.  All  of  the  experimental  studies  cited  here 
agree  well  on  the  value  of  C44,  and  using  EDIP  to  generate 
and  calculate  the  elastic  constants  of  amorphous  samples 
yields  a  result  consistently  about  a  factor  of  2  below  this 
value.  This  may  be  related  to  EDIP’s  neglect  of  rehydridiza- 
tion  for  strains  small  enough  not  to  change  the  coordination, 
which  is  also  responsible  for  EDIP’s  underestimation  of  the 
Kleinman  parameter  and  the  relaxed  C44  of  the  crystal.12  For 
small  distortions  of  the  diamond  lattice  (Z=  4),  EDIP  and  SW 
are  both  equivalent  to  the  rigid  hybrid  approximation,  which 
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TABLE  V.  Elastic  constants  for  the  crystal,  amorphous  sample 
G,  and  the  amorphous  pocket  composite  sample.  The  Young’s 
modulus  value  ( E )  given  for  the  crystal  and  the  amorphous  pocket 
sample  are  the  average  of  the  Reuss  and  Voigt  spatial  averages  (Ref. 
42).  Upper  and  lower  bounds  based  on  Eq.  (9)  are  shown  for 
comparison. 


c„ 

C\2 

C44 

E 

Crystal 

172 

64.7 

72.8 

159 

Amorphous  G 

131 

81.4 

24.8 

68.6 

Composite 

169 

65.4 

69.9 

154 

Upper  bound 

156 

Lower  bound 

151 

is  very  accurate  for  unrelaxed  crystal  elasticity,  but  breaks 
down  with  internal  relaxation.27 

C.  Amorphous  pocket  sample 

The  elastic  constants  of  our  amorphous  pocket  sample  are 
shown  (along  with  those  of  the  component  parts)  in  Table  V. 
It  can  be  seen  that  the  Young’s  modulus  of  the  composite  is 
within  the  upper  and  lower  bounds  predicted  by  Eq.  (9).  This 
was  true,  in  fact,  for  all  of  several  different  methods  we  used 
for  applying  Eq.  (9),  including  using  the  Young’s  modulus 
for  (100)  and  (110)  directions  (as  opposed  to  using  spatially 
averaged  values).  This  result  agrees  with  experimental  obser¬ 
vations  of  the  elastic  properties  of  silicon  during 
amorphization.4 

Predictions  based  on  Eq.  (9)  did  not  hold,  however,  for 
elastic  constants  calculated  using  a  similarly  prepared  com¬ 
posite  sample  that  skipped  the  annealing  step.  The 
crystalline-amorphous  interface,  when  only  relaxed  at  0  K, 
seems  the  likely  reason  for  this. 

D.  Isolated  defects 

The  interstitial  defects  that  resulted  from  the  random- 
placement  generation  process  described  in  Sec.  II  A  1  are 
illustrated  in  Fig.  3.  It  can  be  seen  that  the  defects  DBb  DB2, 
and  DB3  are  minor  variations  on  the  (110)  dumbell  (DB) 
defect.  In  each  case,  the  atom  positions  are  very  similar,  but 
sufficiently  different  to  cause  the  bonding  geometry  to 
change.  All  of  these  configurations  relaxed  to  the  (110)  dum¬ 
bell  configuration  upon  annealing  the  isolated  defect  at  room 
temperature.  We  infer  that  these  defects  were  stabilized  by 
the  presence  of  other  defects. 

The  DB4  defect  complex  is  also  a  variation  on  the  (110) 
dumbell  in  which  two  neighbors  of  the  atom  that  is  “split”  to 
form  the  dumbell  are  moved  together  to  bond.  While  the 
other  dumbell  variations  may  simply  be  EDIP  artifacts,  this 
defect  has  been  seen  before  with  the  Tersoff  3  potential.50 

The  formation  energy  of  3.38  eV  calculated  for  the  (110) 
dumbell  defect  does  not  precisely  agree  with  one  (3.35  eV) 
previously  published.13  This  is  perhaps  partly  due  to  our  in¬ 
dependent  relaxation  of  all  three  cell  dimensions  in  addition 
to  atomic  coordinates,  as  well  as  to  size  effects  with  periodic 
boundary  conditions,  as  suggested  by  Fig.  4.  Effects  of  an 


overlapping  strain  field  due  to  periodic  boundary  conditions 
don’t  appear  to  be  of  concern  (at  least  to  the  level  of  preci¬ 
sion  discussed  here)  for  our  1728  atom  samples.  While  we 
did  not  verify  this  for  each  of  the  many  samples  discussed  in 
the  present  work,  we  take  confidence  in  the  leveling  off  of 
the  curve  in  Fig.  4. 

We  calculated  the  unrelaxed  formation  energies  of  the  te¬ 
tragonal  and  hexagonal  interstitials  for  EDIP,  and  obtained 
10.58  and  6.85  eV,  respectively,  in  agreement  with  Justo  el 
al.  We  noted,  however,  that  these  two  interstitial  positions 
are  not  stable  in  EDIP,  as  shown  in  Fig.  5.  An  interstitial  will 
shift  away  from  these  positions  upon  relaxation.  In  fact,  even 
at  0  K,  numerical  error  generated  during  the  relaxation  of 
atomic  positions  provided  sufficient  asymmetry  to  cause  the 
hexagonal  interstitial  to  relax  into  an  off-center  position 
(hexA).  This  asymmetric  position  remains  sixfold  coordi¬ 
nated.  To  obtain  a  relaxed  formation  energy  for  the  hexago¬ 
nal  interstitial,  therefore,  the  interstitial  position  had  to  be 
artificially  maintained.  This  relaxed  configuration  for  the 
hexagonal  defect  was  metastable,  with  a  formation  energy  of 
4.19  eV.  After  room-temperature  annealing,  both  hexagonal 
and  tetragonal  interstitials  relaxed  into  a  (110)  dumbell.  This 
instability  of  the  tetrahedral  interstitial  agrees  with  LDA 
predictions.51  The  result  that  the  unrelaxed  hexagonal  inter¬ 
stitial  is  unstable,  but  that  relaxing  the  system  while  con¬ 
straining  symmetry  among  its  neighbors  yields  a  metastable 
defect,  also  agrees  with  the  ab  initio  result.52 

We  also  note  in  passing  that  a  bond  defect  involving  no 
coordination  defects  is  stable  using  EDIP,  the  relaxed  forma- 


FIG.  3.  Interstitials  generated  by  random  placement  and  subse¬ 
quent  relaxation  using  EDIP.  The  corners  of  a  diamond  unit  cell  are 
connected  by  lines,  and  atoms  occupying  interstitial  positions  are 
darkened.  Bonds  are  drawn  between  atoms  closer  than  2.56  A.  H2 
involves  two  atoms  which  are  extra  to  the  diamond  unit  cell — all 
other  defects  pictured  here  involve  only  one  extra  atom. 
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FIG.  4.  Formation  energy  for  (110)  dumbell  calculated  in  vari¬ 
ous  sized  supercells  with  periodic  boundary  conditions. 

tion  energy  and  formation  volume  of  which  are  2.38  eV  and 
+  1.1  A3,  respectively.  This  defect  has  been  studied  numerous 
times  under  various  names.  We  adopt  the  latest  nomenclature 
for  this  defect:  the  fourfold  coordinated  defect  (FFCD).53 
The  (FFCD)  is  pictured  in  Fig.  6. 

The  earliest  (to  our  knowledge)  descriptions  of  this  defect 
were  by  Stillinger  and  Weber,9  and  by  Wooten,  Winer,  and 
Weaire,54  who  used  it  to  generate  samples  of  amorphous  Si. 
Later,  Motooka55  described  how  a  divacancy/di-interstitial 
complex  results  in  an  equivalent  configuration.  Tang  et  al.56 
later  rediscovered  this  arrangement  as  the  result  of  the  close 
approach  of  a  vacancy-interstitial  pair,  and  made  a  tight- 
binding  calculation  of  the  formation  energy  (3.51  eV).  Stock 
et  al.51  identified  the  Stillinger  and  Weber  bond  defect,  the 
WWW  bond-switching  mechanism,  and  the  Tang  “I-V  com¬ 
plex”  as  one  and  the  same.  Cargoni  et  al.5&  examined  the 
bond  defect  with  tight-binding  molecular  dynamics  (TBMD) 
and  ab  initio  Hartree-Fock  calculations,  reporting  a  forma¬ 
tion  energy  of  3.26  eV.  Marques  et  al.50  used  the  Tersoff  3 
(T3)  potential11  to  calculate  a  formation  energy  of  3.01  eV. 
Recently,  Goedecker  et  al.53  used  density  functional  theory 
(DFT)  to  calculate  the  formation  energy  of  this  defect.  The 
local  density  approximation  (LDA)  predicted  2.34  eV,  and 
the  general  gradient  approximation  (GGA)  predicted 
2.42  eV.  These  latest  ab  initio  calculations  are  in  remarkable 
agreement  with  EDIR 

It  is  interesting  to  note  that  while  the  formation  energy  of 
this  defect  is  dramatically  lower  than  that  of  the  other  defects 
present  in  our  samples,  it  does  not  occur  in  any  of  them.  This 


FIG.  5.  Unrelaxed  formation  energy  for  an  interstitial  placed 
along  the  unit  cell  body  diagonal  in  Si,  calculated  using  EDIP.  The 
solid  and  dashed  vertical  lines  indicate  the  positions  of  the  hexago¬ 
nal  and  tetragonal  interstials,  respectively. 


FIG.  6.  The  fourfold  coordinated  defect  referred  to  a  unit  cell. 
The  two  displaced  atoms  are  darkened. 

is  perhaps  due  to  the  orchestrated  movement  of  atoms  that  is 
required  to  produce  it,  which  decreases  the  likelihood  of  it 
naturally  occuring. 

Three  types  of  divacancies  cropped  up  in  the  samples 
with  random  vacancy  placement.  One  involved  missing  near¬ 
est  neighbors  {V2N)  with  symmetric  outward  breathing.  An¬ 
other  involved  missing  next-nearest  neighbors  (V^aw)-  The 
third  also  consisted  of  missing  next-nearest  neighbors,  but 
differed  in  that  the  common  neighbor  to  the  two  vacant  sites 
shifted  to  become  threefold  coordinated  ( V2a )■  When  isolated 
in  a  supercell,  V2a  survived  annealing,  but  V2NN  relaxed  into 
a  fourth  configuration,  V^aw n  involving  missing  nearest 
neighbors  that  differed  from  V2N  in  that  instead  of  involving 
six  threefold  coordinated  atoms,  it  had  only  three  threefold 
coordinated  atoms  accompanied  by  a  single  fivefold  coordi¬ 
nated  atom.  These  vacancy  configurations  are  illustrated  in 
Fig.  7.  Though  one  case  of  the  divacancy  V2NN  had  survived 
the  annealing  step  in  one  of  the  random-placement  samples, 
it  was  not  in  isolation  in  the  supercell,  and  the  strain  fields  of 
the  other  defects  present  may  have  served  to  stabilize  it  to 
some  degree. 

The  relaxed  monovacancy  and  divacancy  energies  show 
good  qualitative  agreement  with  DFT/LDA  calculations, 
which  predict  3.29,  4.63,  and  5.90  eV  for  V,  V2N,  and  V2NN, 
respectively.59  The  binding  energy  for  the  divacancy  is  there¬ 
fore  1.60  and  1.95  eV  according  to  EDIP  and  DFT/LDA, 
respectively.  We  note  that  while  the  monovacancy,  “simple” 
divacancy,  and  “split”  divacancy  of  Seong  et  al.59  have  the 
same  unrelaxed  structures  as  V,  V2N,  and  V2NN,  respectively, 
the  energies  we  cite  here  are  for  the  relaxed  structures,  on 
which  DFT/LDA  and  EDIP  differ.  In  the  case  of  V,  for  ex¬ 
ample,  Seong  et  al.  report  a  monovacancy  structure  that  is 
somewhere  between  the  symmetric  V  favored  by  EDIP  and 
the  distorted  VjX. 

The  formation  volume  for  an  isolated  vacancy  was  found 
to  be  +28.8  A3.  The  EDIP  monovacancy  is  therefore 
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FIG.  7.  Divacancy  complexes  referred  to  a  unit  cell  (left  to 
right):  V2aw>  U2 a,  and  V2NNr.  Highlighted  and  darkened  atoms  are 
twofold  and  threefold  coordinated,  respectively.  At  left,  two  next- 
nearest  neighbors  are  missing  from  the  upper  half  of  the  unit  cell. 
At  center,  the  twofold  coordinated  atom  (marked  with  a  cross)  has 
broken  symmetry  and  become  threefold  coordinated.  Annealing 
Vinn  at  room  temperature  yielded  the  rightmost  figure,  in  which  the 
uppermost  atom  (black)  is  fivefold  coordinated.  This  atom  is  not 
part  of  the  unit  cell  shown,  but  is  included  to  provide  context.  While 
V2NNr  was  generated  by  annealing  V^iW"  4  actually  involved  miss¬ 
ing  nearest  neighbors,  similar  to  the  simpler  V2n  (n°t  shown). 

outward-breathing  (the  atomic  volume  of  the  EDIP  crystal  is 
20.0  A3,  in  agreement  with  early  ab  initio  calculations.60,61 
Later  calculations  have  predicted  an  inward-breathing 
vacancy,59,62,63  and  a  recent  DFT/LDA  calculation  study 
showed  that  the  relaxation  around  a  neutral  monovacancy 
has  a  formation  volume  of  -1.7  A3.64  The  negative  forma¬ 
tion  volume  accompanied  a  Jahn-Teller  distortion  in  which 
the  neighbors  of  the  vacant  site  bond  pairwise  across  (110) 
directions.  While  we  found  that  that  such  a  defect  is  stable 
using  EDIP  (and  it  did  indeed  have  a  negative  formation 
volume  of  -14.5  A3),  its  formation  energy  was  higher  than 
that  of  the  monovacancy:  4.06  eV.  When  comparing  forma¬ 
tion  volumes  of  various  defects,  it  is  helpful  to  remember 
that  the  unrelaxed  vacancy  has  a  formation  volume  of 
+20  A3,  whereas  that  of  an  unrelaxed  monointerstitial  con¬ 
figuration  is  -20  A3.  A  negative  formation  volume,  there¬ 
fore,  is  not  as  surprising  for  an  interstitial  as  for  a  vacancy, 
where  relaxation  must  be  very  significant  to  pack  atoms 
more  efficiently  than  the  crystal. 

The  defect  (VJT)  disappeared  upon  annealing  at  300  K  for 
500  ps  in  three  samples  we  tested,  evolving  each  time  into  a 
defect  (VJTr)  involving  one  three-coordinated  and  one  five- 
coordinated  atom  and  having  a  formation  energy  and  volume 
of  3  eV  and  -6.7  A3,  respectively.  We  calculated  the  effect 
of  both  V3j  and  VyYr  on  the  elastic  constants  at  0  K.  Whereas 
the  monovacancy  V  is  seen  to  cause  a  reduction  in  every 
elastic  constant,  the  other  two  monovacancy  types  show  be¬ 
havior  similar  to  the  FFCD  in  that  they  actually  stiffen  Cl2 
and  the  bulk  modulus.  In  each  case,  however,  the  Young’s 
modulus  is  observed  to  lessen.  The  three  monovacancy  con¬ 
figurations  are  illustrated  in  Fig.  8.  The  comparison  of  VjTr 
with  V2NNr  reveals  a  common  structure:  each  has  a  fivefold 
coordinated  atom  with  four  bonds  almost  being  coplanar,  and 
the  fifth  bond  pointing  in  the  general  direction  of  the  three¬ 
fold  coordinated  atoms  (this  is  not  particularly  evident  in 
Figs.  7  and  8 — only  two  of  the  almost  coplanar  bonds  are 
shown  in  each  case). 

The  results  for  formation  energy,  formation  volume,  and 
effects  on  elastic  constants  of  the  isolated  defect  samples  are 
summarized  in  Table  VI.  While  the  sign  on  the  deviation 


FIG.  8.  Stable  EDIP  monovacancy  configurations.  At  left  is  the 
outward-breathing  vacancy  (V)  observed  in  all  randomly  generated 
samples.  Vacancy  neighbors  are  darkened.  The  Jahn-Teller  distor¬ 
tion  with  bonding  between  pairs  of  vacancy  neighbors  is  shown  at 
center  (VjT),  and  the  result  upon  a  room  temperature  annealing  of 
this  defect  is  shown  at  right  ( VjTr ).  The  cross  marks  an  atom  whose 
position  changes  significantly  between  the  three  configurations 
shown,  as  does  the  circle.  V]Tr  contains  one  five-coordinated  atom 
(circle)  and  non  three-coordinated  atom  (back  lower  left  comer). 

from  crystalline  elastic  constants  varied  between  vacancies 
and  interstitials  in  some  cases  (C12  and  B),  all  defects  were 
observed  to  reduce  Young’s  modulus.  When  taken  on  a  per- 
atom  basis,  the  effect  on  the  (100)  Young’s  modulus  is  quite 
similar  among  all  defect  types  presented  in  Table  VI,  being 
bracketed  between  -0.25  and  -0.45%  (aside  from  the  tetra¬ 
hedral  interstitial,  which  is  locally  unstable  according  to  both 
EDIP  and  LDA  calculations). 

While  one  might  expect  that  greater  formation  volume 
magnitude  would  be  accompanied  by  a  larger  formation  en¬ 
ergy  (due  to  larger  strain  fields),  we  found  that  the  relation¬ 
ship  between  formation  volume  and  and  formation  energy 
among  defects  was  not  completely  rigid.  The  two  defects 
with  the  lowest  formation  energies,  H2  and  the  FFCD,  did 
indeed  have  low  formation  volume  magnitudes.  The  relation¬ 
ship  between  V2N  and  V2NNn  however,  is  the  opposite  of 
what  one  would  expect,  as  is  the  large  formation  volume  that 
accompanies  the  most  stable  monovacancy  V.  One  may  also 
wonder  if  the  overall  positive  formation  volume  may  not  be 
driving  the  decrease  in  elastic  modulus.  This  is  clearly  not 
the  case,  as  can  be  seen  by  considering  Table  VI.  For  every 
type  of  point  defect,  having  positive  or  negative  formation 
volume,  the  change  in  Young’s  modulus  is  negative.  The 
changes  in  other  elastic  constants  are  also  uncorrelated  with 
formation  volume.  This  is  to  be  expected,  since  elastic  con¬ 
stants  are  not  strictly  a  function  of  volume,  but  also  of  bond¬ 
ing  topology  and  bond  strength. 

We  examined  the  changes  in  ring  structure  that  accompa¬ 
nied  the  defects  of  Table  VI  in  an  attempt  to  relate  ring 
statistics  and  coordination  numbers  to  the  observed  changes 
in  elastic  constants.  The  interstitial  defects  which  involved 
extra  atoms  did  not  present  any  obvious  pattern  in  this  re¬ 
gard,  but  the  vacancy  defects  did.  Table  VII  shows  the  ring 
and  coordination  figures  for  the  vacancy  and  FFCD  defects. 
A  ring  is  defined  here  as  a  closed  path  which  is  a  series  of 
sequentially  bonded  atoms  without  overlap,  and  a  primitive 
ring  is  a  ring  which  cannot  be  decomposed  into  two  smaller 
rings.65,66  We  determined  primitive  ring  statistics  using  a  re¬ 
cent  effective  ring  search  algorithm.65 

It  can  be  seen  that  with  increased  numbers  of  undercoor¬ 
dinated  atoms,  the  bulk  modulus  is  lessened.  Concurrent 
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TABLE  VI.  Formation  energies,  volumes,  and  percent  change  in  elastic  constants  for  a  single  defect  in  a 
1728-atom  supercell,  calculated  using  EDIP.  The  numbers  in  parenthesis  represent  how  many  of  the  90 
“extra”  atoms  randomly  inserted  into  our  samples  appeared  in  a  particular  interstitial  type.  The  number  of 
vacant  sites  appearing  in  a  particular  vacancy  complex  is  also  shown  in  parenthesis,  also  out  of  a  total  of  90. 


E/ 

EDIP 

,  (eV) 

LDA 

GGA 

Vfi 

EDIP 

A3 

LDA 

c„ 

C*12 

C44 

B 

E{m 

£<no> 

DB 

3.38  (67) 

2.88“ 

3.31“ 

8.7 

-2.0b 

-0.02 

0.72 

-0.32 

0.29 

-0.35 

-0.25 

DB! 

3.55  (4) 

7.0 

-0.04 

0.64 

-0.35 

0.25 

-0.28 

-0.27 

db2 

3.55  (4) 

6.7 

-0.10 

0.59 

-0.38 

0.20 

-0.41 

-0.31 

db3 

3.53  (6) 

6.2 

-0.03 

0.54 

-0.40 

0.22 

-0.35 

-0.28 

db4 

3.50  (1) 

-0.6 

-0.03 

0.56 

-0.36 

0.22 

-0.29 

-0.26 

Tet 

4.10  (0) 

3.43c 

4.07c 

19.6 

0.12 

0.81 

-0.43 

0.41 

-0.18 

-0.24 

Hex 

4.19  (0) 

2.87“ 

3.31“ 

8.3 

-0.11 

0.50 

-0.45 

0.16 

-0.38 

-0.35 

1  lex4 

3.95  (0) 

7.4 

-0.10 

0.63 

-0.50 

0.21 

-0.42 

-0.38 

FFCD 

2.38  (0) 

2.34“ 

2.42“ 

1.1 

5.9b 

-0.12 

0.18 

-0.31 

0.01 

-0.26 

-0.25 

h2 

4.81  (8) 

-4.0 

-0.09 

1.00 

-0.42 

0.43 

-0.58 

-0.36 

V 

3.22  (84) 

3.29d 

28.8 

12. 8d 

-0.31 

-0.42 

-0.27 

-0.34 

-0.26 

-0.28 

PjT 

4.06  (0) 

3.49e 

-14.5 

-1.7e 

-0.21 

0.32 

-0.44 

0.03 

-0.45 

-0.38 

VjTr 

3.65  (0) 

-6.7 

-0.18 

0.24 

-0.32 

0.01 

-0.37 

-0.29 

V2N 

4.84  (2) 

4.63d 

53.4 

16.3d 

-0.54 

-0.61 

-0.46 

-0.55 

-0.50 

-0.49 

V2NN 

6.77  (2) 

5.90d 

51.3 

14.1d 

-0.67 

-0.81 

-0.57 

-0.72 

-0.61 

-0.60 

V2NNr 

5.23  (0) 

16.2 

-0.38 

0.05 

-0.50 

-0.18 

-0.57 

-0.48 

Via 

5.84  (2) 

47.6 

-0.61 

-0.75 

-0.58 

-0.65 

-0.54 

-0.58 

“From  Ref.  53. 
bFrom  Ref.  56. 

“From  Ref.  51. 
dFrom  Ref.  59. 

“Front  Ref.  64. 

with  this  trend  are  the  softening  effects  of  the  loss  of  six- 
rings  and  the  gain  of  larger  rings.  Offsetting  factors  include 
the  gain  of  overcoordinated  atoms  and  rings  of  size  less  than 
6.  The  relationship  between  coordination,  ring  size  popula¬ 
tion,  and  bulk  modulus  was  more  complicated  for  many  of 
the  defects  not  shown  in  Table  VII,  but  the  greatest  stiffening 
of  the  bulk  modulus  did  coincide  with  the  H2  and  tet  defects, 
which  were  characterized  by  overcoordinated  atoms  and 
small  rings.  The  numbers  for  these  other  defects  are  shown 
in  Table  VIII.  The  hex  and  hex4  defects  would  appear  to  be 
likely  to  raise  the  bulk  modulus  more  than  other  defects 
above  them  in  Table  VIII,  given  their  three-rings  and  their 
lack  of  large  rings,  but  perhaps  their  peculiar  planar  geom¬ 
etry  partially  nullifies  these  stiffening  characteristics. 

It  it  interesting  to  note  that  while  EDIP  may  have  erred  in 
predicting  hexA  to  have  a  lower  formation  energy  than  hex, 
the  effects  on  the  supercell  elastic  constants  are  much  the 
same  in  either  case.  The  ordering  of  the  formation  energies 
of  the  tetrahedral  and  the  hex4  interstitials  agrees  with  LDA 
calculations,51  while  this  was  not  the  case  with  the  symmet¬ 
ric  hexagonal  interstitial. 

Finally,  we  used  the  numbers  in  Table  VI  to  predict  the 
results  for  the  samples  (summarized  in  Fig.  1)  generated  by 
random  defect  placement.  We  compared  the  elastic  constant 
changes  in  these  samples  to  the  sum  of  the  individual  con¬ 


tributions  listed  in  Table  VI  for  the  defects  involved.  This 
sum  usually  overestimated  the  change  in  elastic  constants, 
making  the  worst  approximation  for  the  samples  containing 
the  highest  defect  concentration.  The  error  was  in  some  cases 
as  high  as  33%,  but  averaged  only  10%.  The  accumulation  of 
elastic  constant  changes  was  therefore  nearly  linear  with  in¬ 
creasing  defect  content,  as  was  shown  in  Fig.  1,  though  some 
saturation  effects  are  understandably  present  at  such  high 
defect  concentrations. 


IV.  DISCUSSION 
A.  Defects  and  elasticity 

We  have  used  a  well  tested,  physically  motivated  potential 
to  study  the  effects  of  defects  on  the  elastic  constants  of 
silicon.  That  our  results  contradict  those  of  Clark  and 
Ackland32  is  not  surprising,  in  view  of  the  untested  potential 
they  used  and  the  lack  of  essential  physics,  such  as  explicit 
angular  dependence.  The  elastic  constants  of  the  CA  poten¬ 
tial  are  significantly  different  from  experiment,  even  more 
than  the  authors  may  have  realized.  They  appear  to  have 
mistakenly  used  the  experimental  constants  of  Si  at  1477  K 
as  a  reference  for  their  calculations  at  0  K.32,37  The  elastic 
constants  of  the  CA  potential  differ  from  experiment38  by  11, 
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TABLE  VII.  Coordination  and  ring  statistics  changes  compared 
to  the  percent  change  in  bulk  modulus  for  vacancy  defects  and  the 
FFCD.  Net  changes  are  for  isolated  defects  in  a  1728-atom  super¬ 
cell,  referred  to  a  similar  crystalline  supercell.  The  net  change  in  the 
number  of  atoms  having  a  particular  coordination  is  given  by  N. 
The  numbers  under  the  “rings”  heading  refer  to  the  net  gain  or  loss 
of  rings  of  a  particular  size. 


Defect 

N 

Rings 

A B,  % 

2 

3 

5 

5 

6 

7 

9 

11 

12 

14 

Eft 

0 

0 

0 

4 

-12 

0 

4 

0 

0 

0 

0.03 

FFCD 

0 

0 

0 

4 

-12 

8 

0 

0 

0 

0 

0.01 

VjTr 

0 

1 

1 

4 

-10 

0 

2 

0 

0 

0 

0.01 

0 

3 

1 

4 

-16 

0 

0 

2 

0 

0 

-0.18 

V 

0 

4 

0 

0 

-12 

0 

0 

0 

4 

0 

-0.34 

v2.  n 

0 

6 

0 

0 

-18 

0 

0 

0 

2 

9 

-0.55 

v2a 

0 

6 

0 

2 

-22 

0 

0 

0 

7 

0 

-0.65 

V2NN 

1 

6 

0 

0 

-22 

0 

0 

0 

8 

0 

-0.72 

91,  and  67  %,  for  Cu,  C12,  and  C44,  respectively.  The  dis¬ 
agreement  in  results  is  perhaps  best  summarized  by  consid¬ 
ering  Young’s  modulus,  which  is  commonly  of  importance  in 
practical  situations.  We  observed  that  all  defect  types, 
whether  isolated  or  in  random  arrangements,  caused  a  less¬ 
ening  of  Young’s  modulus,  whereas  the  results  of  Clark  and 
Ackland  imply  that  the  presence  of  vacancies  or  interstitials 
causes  Young’s  modulus  to  increase. 

The  only  other  study  of  the  effects  of  defects  on  elastic 
constants  that  we  were  able  to  find  in  the  literature  was  the 
tight-binding  molecular  dynamics  work  of  De  Sandre  et  a/.45 
It  is  difficult  to  compare  our  findings  with  regard  to  point 
defects  with  those  of  that  study  given  that  the  lowest  concen¬ 
tration  considered  there  was  9.3%,  at  which  point  the  crystal- 
to-amorphous  transition  has  begun.67  Nevertheless,  the  au¬ 
thors  report  “the  evolution  of  the  elastic  constants  towards  an 
overall  softening  during  the  crystal-to-amorphous  transi¬ 
tion.” 

Our  results  also  qualitatively  agree  with  the  experimental 
data  available  in  the  literature.  Burnett  and  Briggs  made 
measurements  of  the  elastic  constants  of  silicon  which  had 
been  bombarded  with  arsenic  and  silicon  ions.1  Their 
samples  varied  from  being  near  the  amorphization  threshold 
to  being  completely  amorphized,  but  in  every  case,  both  Cn 
and  C44  decreased  as  the  radiation  damage  increased. 

As  an  additional  check  on  our  results,  we  created  more 
defect-containing  simulation  cells  using  the  Stillinger- Weber 
(SW)  potential9  as  we  had  done  with  EDIP  in  Sec.  IIA1. 
One  set  of  cells  contained  vacancies,  the  other  interstitials, 
both  in  varying  concentrations  as  in  Fig.  1.  In  each  case,  the 
Young’s  modulus  along  (110)  and  (100)  decreased  with  in¬ 
creasing  defect  content.  Furthermore,  an  SW  simulation  of 
the  isolated  defects  V,  DB,  and  tet,  as  in  Table  VI,  predicted 
a  decrease  in  both  En 00)  and  E/u 0)  in  all  three  cases.  Cor¬ 
roboration  by  the  most  widely  used  potential  for  silicon  sug¬ 
gests  that  our  general  conclusion  is  robust. 

Since  it  would  appear  from  Table  III  that,  even  at  the  high 
defect  concentrations  studied  here,  point  defects  have  an  ef- 


PHYSICAL  REVIEW  B  70,  134113  (2004) 


TABLE  VIII.  Coordination  and  ring  statistics  changes  compared 
to  the  percent  change  in  bulk  modulus  for  interstitial  defects,  as  in 
Table  VII. 


Defect 

N 

Rings 

AS,  % 

5 

6 

3 

4 

5 

6 

7 

h2 

8 

0 

6 

0 

2 

-2 

0 

0.43 

Tet 

4 

0 

0 

6 

0 

0 

0 

0.41 

DB 

2 

0 

2 

0 

0 

0 

2 

0.29 

DB! 

2 

0 

2 

0 

1 

-4 

4 

0.25 

DB4 

2 

0 

0 

1 

8 

-12 

8 

0.22 

db3 

4 

0 

3 

1 

0 

-2 

3 

0.22 

I  lex<i 

6 

1 

6 

0 

0 

5 

0 

0.21 

DB2 

2 

0 

2 

0 

1 

-4 

3 

0.20 

Hex 

6 

1 

6 

0 

0 

5 

0 

0.16 

feet  on  the  elastic  constants  that  is  fairly  independent  of  one 
another,  the  details  of  coordination  within  each  defect  as  it 
relates  to  connectivity  in  the  lattice  appears  to  be  important. 
As  shown  in  Table  VI,  a  vacancy  will  have  different  effects 
on  the  elastic  constants  depending  on  its  relaxed  bonding 
configuration.  While  EDIP  does  (by  design)  produce  elastic 
constants  for  silicon  close  to  those  of  experiment,  the  con¬ 
figurations  of  relaxed  point  defects  it  produces  do  not  corre¬ 
spond  precisely  with  those  predicted  by  ab  initio  calcula¬ 
tions. 

We  should  mention  that  the  volume  expansion  shown  in 
Fig.  2  for  EDIP  is  likely  not  in  good  agreement  with  reality, 
since  the  samples  involved  contained  mostly  monovacancies, 
which  were  shown  here  to  be  outward  breathing  in  EDIP.  As 
noted  above,  recent  ab  initio  calculations  predict  the  the 
monovacancy  to  have  a  negative  formation  volume,  which 
would  correspond  to  a  densification  with  increasing  monova¬ 
cancy  content.  We  must  therefore  be  cautious  as  to  which 
results  presented  here  we  expect  to  correspond  to  reality. 

This  word  of  caution  does  not  preclude  two  generaliza¬ 
tions  of  practical  interest,  however.  We  reiterate  that  the  gen¬ 
eral  trend  of  Young’s  modulus  is  consistently  downward  with 
increasing  defect  content,  regardless  of  the  positive  or  nega¬ 
tive  formation  volume  of  the  defects  involved.  This  effect 
varies  within  a  very  narrow  range  for  the  defect  configura¬ 
tions  considered  here.  It  should  therefore  be  possible  to  make 
reasonable  predictions  of  the  shift  in  the  Young’s  modulus  of 
a  radiation-damaged  silicon  sample  in  which  isolated  point 
defects  dominate,  even  if  the  exact  geometry  of  each  point 
defect  is  not  well  known.  The  calculation  presented  here  of 
the  crystalline/amorphous  composite  sample  is  also  encour¬ 
aging.  It  appears  that  a  reasonable  prediction  of  changes  in 
Young’s  modulus  can  also  be  made  for  radiation-damaged 
samples  in  which  amorphous  regions  are  important,  based 
simply  on  the  volume  fraction  of  amorphous  material,  ac¬ 
cording  to  Eq.  (9). 

B.  Interatomic  potentials 

There  are  advantages  to  using  an  empirical  potential,  in¬ 
stead  of  ab  initio  methods,  for  the  work  presented  here.  The 


134113-11 


ALLRED  et  al. 


PHYSICAL  REVIEW  B  70,  134113  (2004) 


first  is  that  since  empirical  potentials  will  always  be  compu¬ 
tationally  cheaper  than  ab  initio  calculations,  it  is  important 
to  continue  the  development  and  characterization  of  such 
potentials  in  order  to  enable  the  study  of  systems  many  times 
the  size  of  those  accessible  to  ab  initio  methods.  Secondly, 
despite  the  fact  that  powerful  computers  could  perform  ab 
initio  calculations  on  the  1728  atom  cells  (the  largest  we 
considered),  the  large  number  of  such  cells  and  the  demands 
of  the  numerous  annealing  steps  involved  in  the  preparation 
of  our  simulation  cells  would  have  been  impractical.  We  pre¬ 
pared  over  60  cells  containing  1728  atoms,  each  one  requir¬ 
ing  more  than  250  K  computational  steps  for  the  preparation 
alone.  The  use  of  an  empirical  potential  allowed  us  to  com¬ 
pile  statistics  on  the  effects  of  random  defect  arrangements 
and  notice  trends,  something  that  would  not  have  been  pos¬ 
sible  using  ab  initio  methods. 

We  briefly  discuss  what  our  results  might  imply  about  the 
physics  of  covalent  bonding.  The  fact  that  we  get  qualita¬ 
tively  similar  results  with  ED  IP  and  SW,  which  differ  sharply 
with  CA,  confirms  the  common  wisdom  that  angular  forces 
are  present  in  covalent  solids,  as  argued  by  Born  almost  a 
century  ago.  The  pairwise  repulsion  of  first  neighbors  in  the 
CA  potential,  similar  to  the  early  potential  of  Pearson  et  al.(M 
[which  also  does  not  perform  as  well  as  SW  (Ref.  29)],  is  not 
equivalent  to  an  explicit  angular  force. 

For  the  reasons  mentioned  in  the  introduction,  however, 
the  reasonable  results  of  SW  for  defects  are  quite  surprising 
because  the  potential  can  only  be  justified  for  weakly  dis¬ 
torted,  rigid  sp 3  hybrids.  There  is  no  question  that  both  the 
bond  order26  and  the  strength  of  angular  forces12  depend  on 
the  local  atomic  environment,  as  also  implied  by  analytical 
approximations  of  tight-binding  models.28  EDIP  attempts  to 
capture  these  effects  with  only  a  simple  scalar,  the  coordina¬ 
tion  number.  Of  course,  this  precludes  the  possibility  of 
gradual  rehybridization  under  shear,  before  changes  in  bond¬ 
ing  topology  occur,  and  it  seems  this  is  related  to  the  under¬ 
estimate  of  C 44  for  both  crystal  and  amorphous  structures. 
The  SW  potential  misses  all  of  this  physics,  and  yet  some¬ 
how  manages  to  describe  defect  configrations  fairly  well. 
The  fortuitous  reason  may  be  that  the  “optimal”  angular 
function  for  silicon,  obtained  by  inversion  of  ab  initio  cohe¬ 
sive  energy  curves  for  silicon,  neglecting  environment  de¬ 
pendence,  happens  to  match  the  SW  angular  term  quite 
well.26 

In  any  case,  it  is  clear  that  any  potential  for  silicon  must 
take  into  account  the  local  atomic  environment,  if  not  explic¬ 
itly,  as  in  EDIP,  then  at  least  implicitly  by  considering  all 
neighboring  atoms,  as  in  SW  and  Tersoff.  The  CA  potential. 


similar  to  the  Ackland  potentials  before  it,34,35  makes  the  ad 
hoc  assumption  of  four  diamond-like  bonds  per  atom,  which 
surely  does  not  apply  to  most  defects.  Ackland  acknowl¬ 
edged  the  need  for  coordination  dependence  in  his  original 
paper,34  by  increasing  the  bond  order  for  undercoordinated 
atoms  (Z=  3),  albeit  to  a  value  inappropriate  for  ab  initio 
graphitic  silicon;26  however,  this  modification  was  not  in¬ 
cluded  in  the  CA  potential.32  The  significant  decrease  in 
bond  order  of  roughly  Z_1/2  for  overcoordination  (Z>4)12,26 
was  also  neglected,  so  the  different  conclusion  regarding  the 
effect  of  point  defects  on  elastic  constants  versus  the  present 
study  can  be  attributed  to  an  inadequate  description  of  cova¬ 
lent  bonding. 

V.  CONCLUSIONS 

We  have  characterized  the  structure  of  various  point  de¬ 
fects  according  to  EDIP  and  calculated  their  individual  ef¬ 
fects  on  elastic  constants.  We  have  shown  that,  in  general, 
the  elastic  constants  of  silicon  vary  with  defect  concentration 
in  a  roughly  linear  fashion  at  defect  concentrations  up  to 
0.3%.  We  have  also  demonstrated  the  effects  of  an  amor¬ 
phous  region  embedded  in  a  crystalline  supercell  to  be  well 
described  by  a  simple  equation  involving  the  volume  fraction 
of  amorphous  material.  We  have  also  calculated  elastic  con¬ 
stants  for  several  amorphous  samples,  and  concluded  that,  in 
spite  of  a  good  overall  description  of  the  amorphous  phase, 
EDIP  consistently  underestimates  C44.  As  in  the  case  of  the 
crystal,  this  may  be  due  to  its  neglect  of  rehybridization  un¬ 
der  strain  (as  in  SW). 

The  defect  structures  and  formation  energies  predicted  by 
the  potential,  such  as  the  energetically  favored  fourfold  co¬ 
ordinated  defect,  seem  fairly  realistic,  so  we  may  place  some 
confidence  in  their  effects  on  elasticity.  We  predict  that  both 
interstitial  and  vacancy  defects,  singly  and  in  random  com¬ 
bination,  lessen  the  Young’s  modulus  of  silicon.  We  also  find 
that  the  individual  contributions  of  various  defect  types  to 
changes  in  Young’s  modulus  are  confined  to  a  surprisingly 
small  range.  These  appear  to  be  robust  conclusions  for  po¬ 
tentials  taking  into  account  explicit  angular  forces  for  cova¬ 
lent  bonds. 
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